This is an example of a working script to simulate the 2 stage hard X-ray self-seeding with '400' diamond reflection tuned to 8980.2 eV.


In [ ]:
'''
hard x-ray self-seeding
launch as python hxrss_max.py <Electron Energy in GeV> <Photon Energy in eV> 
'''


import sys, os
import matplotlib.pyplot as plt
from ocelot.gui.accelerator import *
from ocelot.utils.xfel_utils import *
from ocelot.utils.hxrss_common import *
from copy import copy
from ocelot.utils.sim_info import *
import numpy.fft as fft
from ocelot.gui.genesis_plot import *
import time

#####################################################
###########                            ##############
###########         BEGIN OF           ##############
###########  USER-DEFINED QUANTITIES   ##############
###########                            ##############
#####################################################


'''
Notes:

-Current example: 17.5 9000.0
-different nodes allocation changes, in general, the effective ip-seed. 
-For this example to run use salloc -N 5 -p exfel-theory --constraint=AMD --time=120 (even changing the #nodes changes effectively the ip-seed!)
-Time for run: 4052s; time for pp: 119 s
'''

rep_dir   = '/data/netapp/xfel/gianluca/products/desy/xfel/' #repository with the undulator configuration is stored. In this example the configuration file is /data/netapp/xfel/gianluca/products/desy/xfel/sase1/sase1.py
param_dir = '/data/netapp/xfel/gianluca/beams/' #directory where the beam file is stored
exp_dir = '/data/netapp/xfel/gianluca/20pC_9keV_17.5GeV_new/' #directory where the output will be stored
res_dir = exp_dir + 'results/' #Directory where the postprocessing output will be stored
sys.path.append(rep_dir + 'sase1/') 
from sase1 import *



#beamf = param_dir + 'beam_0.1nC_wake.txt'
beamf = param_dir + 'beam_0.02nC_wake_vitali.txt' #beamfile 

run_ids = xrange(0,2) #statistical runs from 0 to n-1 original (0,2)

debug = True
outfig = False
xt_couple = True #Spatio-temporal coupling flag

start_stage = 1
stop_stage = 5
beta_av=25.0

#First HXRSS filter
E_filt_2 = 8980.2 
delay_filt_2 = -2050#-0.00001035

#Second HXRSS filter
E_filt_4 = 8980.2
delay_filt_4 = -2050 #-0.00001035

#Length of undulator parts
z_stop_1 = 6.05*5.0
z_stop_3 = 6.05*5.0
z_stop_5 = 6.05*6.0

#simulation parameters
ncar =151
npart = 8192
dgrid = 4.0e-4

#self-seeding chicanes (already defined by default in sase1.py)
chicane1 = Drift(l=5.1)
chicane1.cryst = Crystal(r=[0,0,0*cm], size=[5*cm,5*cm,100*mum], no=[0,0,-1])
chicane1.cryst.lattice =  CrystalLattice('C')
chicane1.cryst.psi_n = -pi/2. #input angle psi_n according to Authier (symmetric reflection, Si)
chicane1.cryst.ref_idx = (4,0,0)

chicane2 = Drift(l=5.1)
chicane2.cryst = Crystal(r=[0,0,0*cm], size=[5*cm,5*cm,100*mum], no=[0,0,-1])
chicane2.cryst.lattice =  CrystalLattice('C')
chicane2.cryst.psi_n = -pi/2. #input angle psi_n according to Authier (symmetric reflection, Si)
chicane2.cryst.ref_idx = (4,0,0)


#####################################################
###########                            ##############
###########         END OF             ##############
###########  USER-DEFINED QUANTITIES   ##############
###########                            ##############
#####################################################


t0 = time.time()

launcher = get_genesis_launcher()

create_exp_dir(exp_dir, run_ids)
cmd = 'mkdir '+res_dir
if (os.path.isdir(res_dir)==0): os.system(cmd)    
cmd = 'cp '+os.path.realpath(__file__)+' '+exp_dir+'exec_script_'+os.path.basename(__file__)+' '    # write copy of this script to the exp_dir
print cmd                                                                                          
os.system(cmd)                                                                                      

beamg =read_beam_file(beamf)

beam.E = float(sys.argv[1])   # e.g. 14.0
E_ev   = float(sys.argv[2])   # e.g. 8000.0
zmax, Imax = peaks(beamg.z, beamg.I, n=1)
idx = beamg.z.index(zmax)
beamg.idx_max=idx

#overwrites values in /desy/xfel/sas1/sase1.py with data read from beam file
#Careful: beam is of class Beam() while beamg is of class GenesisBeamDefinition()
beam.emit_xn = beamg.ex[idx]
beam.emit_yn = beamg.ey[idx]
beam.gamma_rel = beam.E / (0.511e-3)
beam.emit_x = beamg.ex[idx] / beam.gamma_rel
beam.emit_y = beamg.ey[idx] / beam.gamma_rel
beam.beta_x = beamg.betax[idx]
beam.beta_y = beamg.betay[idx]
beam.alpha_x = beamg.alphax[idx]
beam.alpha_y = beamg.alphay[idx]
beam.C = 0.1
beam.I = beamg.I[idx]
beam.tpulse = (beamg.z[-1]-beamg.z[0])/(299792458.0/1e15)/8.0 #should be using the right speed of light; see line 1405 in adaptors/genesis.py
en_beamfile = beamg.g0[idx]*0.511e-3
beam.emit = ''



lat = MagneticLattice(sase1_segment(n=20))
rematch(beta_av, l_fodo, qdh, lat, extra_fodo, beam, qf, qd) # jeez...

tw0 = Twiss(beam)
tws=twiss(lat, tw0, nPoints = 100) # to make sure the average beta exists, show twiss if needed 

# calculate UR parameters (required later for input generation)
up = UndulatorParameters(und,beam.E)
up.printParameters()

#tapering functions: no taper


taper_func_1 = lambda n : f1(n, [0, 20], 1.0, [0,0], [0,0])
taper_func_2 = lambda n : f1(n, [0, 20], 1.0, [0,0], [0,0])

n0 = [0,8,29,36]
a0 = 1.0
a1 = [0.0, -0.0006,  -0.0006 ]
a2 = [0.0, -0.00013, -0.0001 ]
taper_func_3 = lambda n : f1(n, n0, a0, a1, a2 )

lat1=lat
lat2=lat
lat3=lat

inp = generate_input(up, beam, itdp=True)
maxslice = inp.nslice

inp.lattice_str = generate_lattice(lat1, unit = up.lw, energy = beam.E )
beam_new = transform_beam_file(beam_file = beamf,transform = [ [beam.beta_x,beam.alpha_x], [beam.beta_y,beam.alpha_y] ], energy_scale = beam.E/en_beamfile, emit_scale = 1.0, n_interp=maxslice)
beam_new.gamma_rel = beam.E / (0.511e-3)

inp.beamfile = 'tmp.beam'
inp.beam_file_str = beam_new.f_str= beam_file_str(beam_new)

# plotting beam not required
if debug: plot_beam(plt.figure(), beam_new)
idx_max = beam_new.idx_max
if debug: print 'idx_max', idx_max
if debug: plt.show()

sim_info = SimInfo()
sim_info.log_dir = exp_dir

#################
# stage 1 (SASE)#
#################

if start_stage <= 1 and stop_stage >= 1:
    for run_id in run_ids:
        inp.runid = run_id
        inp.run_dir = exp_dir + 'run_' + str(inp.runid)
        inp.ipseed = 61*(run_id + 1 )

        inp.ncar = ncar
        inp.zstop = z_stop_1
        inp.npart = npart
        inp.dgrid = dgrid
        inp.delz = 1.0

        #inp.nslice=0
        ######inp.zsep = int(beam_new.zsep / inp.xlamds)
        inp.nslice=0
        inp.nharm = 1
        inp.idmppar = 0  
        inp.idmpfld = 1      
        g = run(inp, launcher, readall=True)
        print(g.sliceKeys)
        os.system('cp '+inp.run_dir+'/tmp.beam'+' '+inp.run_dir+'/beam.in.txt')
        #if debug: show_output(g, show_field = True, show_slice=maxslice)
        #if debug: plt.show()
        #log_info(sim_info, g, run_id, 1)
        beam_stage1 = deepcopy(beam_new)
        update_beam_2(beam_stage1, g, maxslice)
        open(inp.run_dir + '/beam.s1.txt','w').write( beam_file_str(beam_stage1) )
        checkout_run(inp.run_dir, run_id, '', '.s1', True)
        #if debug: plt.plot(beam_stage1.dg)
        #if debug: plt.show()
        cmd1 = 'cp '+inp.run_dir+'/tmp.gen'+' '+inp.run_dir+'/geninp.'+str(run_id)+'.s1.inp'
        cmd2 = 'cp '+inp.run_dir+'/lattice.inp'+' '+inp.run_dir+'/lattice.'+str(run_id)+'.s1.inp'
        os.system(cmd1)
        os.system(cmd2)

        if outfig:

            gen_outplot(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s1.gout',save='png',show=False,all=True)
            data_t = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s1.gout',save='png',show=False,far_field=0, freq_domain=0,auto_zoom=0)
            data_f = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s1.gout',save='png',show=False,far_field=0, freq_domain=1,auto_zoom=0)

####################
# stage 2 (filter) #
####################

#-it must be run a first time in debug mode; 
#-then delay_slice and E_ev can be inserted by user
if start_stage <= 2 and stop_stage >= 2:
    print 'Start STAGE 2...'
    #debug = False
    print beam_new.z[-1]
    print beam_new.z[idx_max]
    
    ds = beam_new.z[1]-beam_new.z[0]


    
    E_ev  = E_filt_2  #7995.5#7979.3 #E_ev # To be inserted; this is now the seed ev
    delay_slice = (delay_filt_2 -idx_max)#-14.23e-6#(-3854 -idx_max) #negative is a delay; subtracting idx_max we obtain the delay to have the wanted point aligned with max curr
    #delay_slice = int((delay_filt_2 -np.abs(beam_new.z[idx_max]-beam_new.z[0]))/ds)
   
    for run_id in run_ids:
        inp.runid = run_id
        inp.run_dir = exp_dir + 'run_' + str(inp.runid)
        output_files = [
            inp.run_dir + '/run.' + str(inp.runid) + '.s2.gout.dfl',
            res_dir+'s2.Tmod.'+str(run_id)+'.dat',
            res_dir+'s2.Tpha.'+str(run_id)+'.dat',
            res_dir+'s2.Tmod_large.'+str(run_id)+'.dat',
            res_dir+'s2.Tpha_large.'+str(run_id)+'.dat',
            #res_dir+'s1.Pout.'+str(run_id)+'.dat',
            #res_dir+'s1.Ppha.'+str(run_id)+'.dat',
            #res_dir+'s1.Sout.'+str(run_id)+'.dat',
            #res_dir+'s1.Spha.'+str(run_id)+'.dat',
            res_dir+'s2.Pout_large.'+str(run_id)+'.dat',
            res_dir+'s2.Ppha_large.'+str(run_id)+'.dat',
            res_dir+'s2.Sout_large.'+str(run_id)+'.dat',
            res_dir+'s2.Spha_large.'+str(run_id)+'.dat',
            #res_dir+'s2.Pout.'+str(run_id)+'.dat',
            #res_dir+'s2.Ppha.'+str(run_id)+'.dat',
            #res_dir+'s2.Sout.'+str(run_id)+'.dat',
            #res_dir+'s2.Spha.'+str(run_id)+'.dat'
            ]            
        run_dir = exp_dir + 'run_' + str(run_id)
        input_file = run_dir + '/run.' + str(run_id) + '.s1.gout'
        ssc, Pout, lsc, Sout = sseed_2(input_file=input_file, output_files=output_files, E_ev=E_ev, chicane = chicane1, run_dir=run_dir, delay = delay_slice, debug=debug, xt_couple = xt_couple)        		
        if run_id == run_ids[0]:
            P_av = Pout
            S_av = Sout
        else:
            P_av = P_av + Pout
            S_av = S_av + Sout
    P_av = P_av/len(run_ids)
    S_av = S_av/len(run_ids)
    out = readGenesisOutput(inp.run_dir+'/run.'+str(inp.runid)+'.s1.gout',readall=0,debug=0)
    nsli = out.nSlices
    nall = len(P_av)
    zscale = np.arange(nall) - (nall/2-nsli/2)

    f1 = open(res_dir+'s2.Pout_large_ave.dat','w')
    f2 = open(res_dir+'s2.Sout_large_ave.dat','w')
    f3 = open(res_dir+'s2.Delay_Pout_ave.dat','w')
    for i in range(len(ssc)):
        f1.write('%s ' %(ssc[i]) + '%s '%(P_av[i])+'\n')
        f2.write('%s ' %(lsc[i]) + '%s '%(S_av[i])+'\n')
        f3.write('%s ' %(zscale[i]) + '%s '%(P_av[i])+'\n')
    f1.close()
    f2.close()
    f3.close()
    
    if debug:        
        plt.figure()    
        plt.yscale('log')    
        plt.plot(ssc, P_av)	
        plt.figure()
        plt.plot(h*c/lsc, S_av)	
        plt.figure()
        plt.yscale('log')
        plt.plot(zscale,P_av)
        plt.show()
        dfl = readRadiationFile(inp.run_dir+'/run.'+str(inp.runid)+'.s2.gout.dfl',out.ncar)
        gen_outplot_dfl(dfl,out,save=0,show=True,far_field=0,auto_zoom=0,xy_lim=[])      
        dfl1 = readRadiationFile(inp.run_dir+'/run.'+str(inp.runid)+'.s1.gout.dfl',out.ncar)
        gen_outplot_dfl(dfl1,out,save=0,show=True,far_field=0,auto_zoom=0,xy_lim=[])


########################
# stage 3 (seed ampli) #
########################
    
# stage 3
if start_stage <= 3 and stop_stage >= 3:
               
    deltagave = 0.0
    beam_new = transform_beam_file(beam_file = beamf,transform = [ [beam.beta_x,beam.alpha_x], [beam.beta_y,beam.alpha_y] ], energy_scale = beam.E/en_beamfile, emit_scale = 1.0, n_interp=maxslice)   
    a = beam_new.g0[beam_new.idx_max]
    for run_id in run_ids:                  
        #note the energy scale: for beamf (previous is divided by 14, beacuse that is the initial value; for the new is 1.0 becaus eis the next
        beam_new = transform_beam_file(beam_file = exp_dir + 'run_' + str(run_id) + '/beam.s1.txt',transform = [ [beam.beta_x,beam.alpha_x], [beam.beta_y,beam.alpha_y] ], energy_scale = 1.0, emit_scale = 1.0, n_interp=maxslice)	    
        beam_new.gamma_rel = beam.E / (0.511e-3)
        b = beam_new.g0[beam_new.idx_max]
        deltagave = deltagave + (b-a)

    deltagave = deltagave/len(run_ids)
    factK = (deltagave*(2+und.Kx*und.Kx)/(und.Kx*a)+und.Kx)/und.Kx
    print 'factK ',factK          
    und.Kx = und.Kx * factK
    up = UndulatorParameters(und)
    up.printParameters()
    inp.lattice_str = generate_lattice(lat2, unit = up.lw, energy = beam.E ) 
        
    
    #debug = False
    for run_id in run_ids:        
        inp.runid = run_id
        inp.run_dir = exp_dir + 'run_' + str(inp.runid)
        inp.fieldfile='run.' + str(inp.runid) + '.s2.gout.dfl'
        inp.ipseed = 17*(run_id + 1 )

        os.system('cp '+inp.run_dir+'/run.' + str(inp.runid) + '.s1.gout '+inp.run_dir+'/run.' + str(inp.runid) + '.gout')

        inp.ncar = ncar
        inp.zstop = z_stop_3
        inp.npart = npart 
        inp.dgrid = dgrid
        inp.delz = 1.0
        inp.nslice=0
        inp.zsep = int(beam_new.zsep / inp.xlamds)
        inp.nharm = 1
        inp.idmppar = 0
        inp.idmpfld = 1 
        inp.beamfile = 'beam.s1.txt'
        inp.beam_file_str = beam_new.f_str = open(inp.run_dir + '/beam.s1.txt').read()       	
        
        # plotting beam not required
        if debug: plot_beam(plt.figure(), beam_new)
        idx_max = beam_new.idx_max
        if debug: print 'idx_max', idx_max
        if debug: plt.show()
        
        g = run(inp, launcher, readall = True)
        #if debug: show_output(g, show_field = True, show_slice=maxslice)
        #if debug: plt.show()
        log_info(sim_info, g, run_id, 3)
        beam_stage3 = deepcopy(beam_new)
        update_beam_2(beam_stage3, g, maxslice)
        open(inp.run_dir + '/beam.s3.txt','w').write( beam_file_str(beam_stage3) )
        checkout_run(inp.run_dir, run_id, '', '.s3', True)
        cmd1 = 'cp '+inp.run_dir+'/tmp.gen'+' '+inp.run_dir+'/geninp.'+str(run_id)+'.s3.inp'
        cmd2 = 'cp '+inp.run_dir+'/lattice.inp'+' '+inp.run_dir+'/lattice.'+str(run_id)+'.s3.inp'
        os.system(cmd1)
        os.system(cmd2)

        if outfig:

            gen_outplot(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s3.gout',save='png',show=False,all=True)
            data_t = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s3.gout',save='png',show=False,far_field=0, freq_domain=0,auto_zoom=0)
            data_f = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s3.gout',save='png',show=False,far_field=0, freq_domain=1,auto_zoom=0)



####################
# stage 4 (filter) #
####################


#-it must be run a first time in debug mode; 
#-then delay_slice and E_ev can be inserted by user
if start_stage <= 4 and stop_stage >= 4:
    print 'Start STAGE 4...'
    
    ds = beam_new.z[1]-beam_new.z[0]
    #debug = False
    E_ev  = E_filt_4 #7995.5 #E_ev # To be inserted; this is now the seed ev
    print beam_new.z[idx_max]-beam_new.z[0]
    delay_slice = (delay_filt_4 -idx_max)
    #delay_slice = int((delay_filt_4 -np.abs(beam_new.z[idx_max]-beam_new.z[0]))/ds)
    
    
    for run_id in run_ids:
        inp.runid = run_id
        inp.run_dir = exp_dir + 'run_' + str(inp.runid)
        output_files = [
            inp.run_dir + '/run.' + str(inp.runid) + '.s4.gout.dfl',
            res_dir+'s4.Tmod.'+str(run_id)+'.dat',
            res_dir+'s4.Tpha.'+str(run_id)+'.dat',
            res_dir+'s4.Tmod_large.'+str(run_id)+'.dat',
            res_dir+'s4.Tpha_large.'+str(run_id)+'.dat',
            #res_dir+'s3.Pout.'+str(run_id)+'.dat',
            #res_dir+'s3.Ppha.'+str(run_id)+'.dat',
            #res_dir+'s3.Sout.'+str(run_id)+'.dat',
            #res_dir+'s3.Spha.'+str(run_id)+'.dat',
            res_dir+'s4.Pout_large.'+str(run_id)+'.dat',
            res_dir+'s4.Ppha_large.'+str(run_id)+'.dat',
            res_dir+'s4.Sout_large.'+str(run_id)+'.dat',
            res_dir+'s4.Spha_large.'+str(run_id)+'.dat',
            #res_dir+'s4.Pout.'+str(run_id)+'.dat',
            #res_dir+'s4.Ppha.'+str(run_id)+'.dat',
            #res_dir+'s4.Sout.'+str(run_id)+'.dat',
            #res_dir+'s4.Spha.'+str(run_id)+'.dat'
            ]            
        run_dir = exp_dir + 'run_' + str(run_id)
        input_file = run_dir + '/run.' + str(run_id) + '.s3.gout'   
        ssc, Pout, lsc, Sout = sseed_2(input_file=input_file, output_files=output_files, E_ev=E_ev, chicane = chicane2, run_dir=run_dir, delay = delay_slice, debug=debug, xt_couple = xt_couple)

        if run_id == run_ids[0]:
            P_av = Pout
            S_av = Sout
        else:
            P_av = P_av + Pout
            S_av = S_av + Sout
    P_av = P_av/len(run_ids)
    S_av = S_av/len(run_ids)
    out = readGenesisOutput(inp.run_dir+'/run.'+str(inp.runid)+'.s3.gout',readall=0,debug=0)
    nsli = out.nSlices
    nall = len(P_av)    
    zscale = np.arange(nall) - (nall/2-nsli/2)
    
    f1 = open(res_dir+'s4.Pout_large_ave.dat','w')
    f2 = open(res_dir+'s4.Sout_large_ave.dat','w')
    f3 = open(res_dir+'s4.Delay_Pout_ave.dat','w')
    for i in range(len(ssc)):
        f1.write('%s ' %(ssc[i]) + '%s '%(P_av[i])+'\n')
        f2.write('%s ' %(lsc[i]) + '%s '%(S_av[i])+'\n')
        f3.write('%s ' %(zscale[i]) + '%s '%(P_av[i])+'\n')
    f1.close()
    f2.close()
    f3.close()
    if debug:
        plt.figure()    
        plt.yscale('log')    
        plt.plot(ssc, P_av)	
        plt.figure()
        plt.plot(h*c/lsc, S_av)
        plt.figure()
        plt.yscale('log')
        plt.plot(zscale,P_av)
        plt.show()
        dfl = readRadiationFile(inp.run_dir+'/run.'+str(inp.runid)+'.s4.gout.dfl',out.ncar)
        gen_outplot_dfl(dfl,out,save=0,show=True,far_field=0,auto_zoom=0,xy_lim=[])	
        dfl1 = readRadiationFile(inp.run_dir+'/run.'+str(inp.runid)+'.s3.gout.dfl',out.ncar)
        gen_outplot_dfl(dfl1,out,save=0,show=True,far_field=0,auto_zoom=0,xy_lim=[])

########################
# stage 5 (seed ampli) #
########################
    

if start_stage <= 5 and stop_stage >= 5:
        
    deltagave = 0.0
    beam_new = transform_beam_file(beam_file = beamf,transform = [ [beam.beta_x,beam.alpha_x], [beam.beta_y,beam.alpha_y] ], energy_new = beam.E, emit_scale = 1.0, n_interp=maxslice)   
    a = beam_new.g0[beam_new.idx_max]
    for run_id in run_ids:           
        beam_new = transform_beam_file(beam_file = exp_dir + 'run_' + str(run_id) + '/beam.s3.txt',transform = [ [beam.beta_x,beam.alpha_x], [beam.beta_y,beam.alpha_y] ], energy_scale = 1.0, emit_scale = 1.0, n_interp=maxslice)	    
        beam_new.gamma_rel = beam.E / (0.511e-3)
        b = beam_new.g0[beam_new.idx_max]
        deltagave = deltagave + (b-a)
    deltagave = deltagave/len(run_ids)
    print deltagave*(2+und.Kx*und.Kx)/(und.Kx*a)
    factK = (deltagave*(2+und.Kx*und.Kx)/(und.Kx*a)+und.Kx)/und.Kx
    print factK          
    und.Kx = und.Kx * factK
    up = UndulatorParameters(und)
    up.printParameters()
    #inp = generate_input(up, beam, itdp=True)
    inp.lattice_str = generate_lattice(lat3, unit = up.lw, energy = beam.E ) 
    ####und.Kx = und.Kx * 1.002 #fine tuning
    
    debug = False
    for run_id in run_ids:        
        inp.runid = run_id
        inp.run_dir = exp_dir + 'run_' + str(inp.runid)
        inp.fieldfile='run.' + str(inp.runid) + '.s4.gout.dfl'
        inp.ipseed = 71*(run_id + 1 )

        os.system('cp '+inp.run_dir+'/run.' + str(inp.runid) + '.s3.gout '+inp.run_dir+'/run.' + str(inp.runid) + '.gout')

        inp.ncar = ncar
        inp.zstop = z_stop_5 
        inp.npart = npart 
        inp.dgrid = dgrid
        inp.delz = 1.0
        inp.nslice=0
        #inp.zsep = int(beam_new.zsep / inp.xlamds)
        inp.nharm = 1
        inp.idmppar = 0   
        inp.idmpfld = 1         
        inp.beamfile = 'beam.s3.txt'
        inp.beam_file_str = beam_new.f_str = open(inp.run_dir + '/beam.s3.txt').read()       

        # plotting beam not required
        if debug: plot_beam(plt.figure(), beam_new)
        idx_max = beam_new.idx_max
        if debug: print 'idx_max', idx_max
        if debug: plt.show()

        g = run(inp, launcher, readall=True)
    
        if debug: show_output(g, show_field = True, show_slice=maxslice)
        if debug: plt.show()
        log_info(sim_info, g, run_id, 5)
        beam_stage5 = deepcopy(beam_new)
        update_beam_2(beam_stage5, g, maxslice)
        open(inp.run_dir + '/beam.s5.txt','w').write( beam_file_str(beam_stage5) )
        checkout_run(inp.run_dir, run_id, '', '.s5', True)
        cmd1 = 'cp '+inp.run_dir+'/tmp.gen'+' '+inp.run_dir+'/geninp.'+str(run_id)+'.s5.inp'
        cmd2 = 'cp '+inp.run_dir+'/lattice.inp'+' '+inp.run_dir+'/lattice.'+str(run_id)+'.s5.inp'
        os.system(cmd1)
        os.system(cmd2)
        
        if outfig:
            gen_outplot(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s5.gout',save='png',show=False,all=True)
            data_t = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s5.gout',save='png',show=False,far_field=0, freq_domain=0,auto_zoom=0)
            data_f = gen_outplot_dfl(exp_dir + 'run_' + str(inp.runid)+'/run.' + str(inp.runid) + '.s5.gout',save='png',show=False,far_field=0, freq_domain=1,auto_zoom=0)

t1 = time.time()
delta = t1-t0
print 'total time = ', delta
print 'Starting post-processing...'
t3 = time.time()
gen_stat_plot(proj_dir=exp_dir,run_inp=[],stage_inp=[],param_inp=[],s_param_inp=['p_int','energy','r_size_weighted'],z_param_inp=['p_int','phi_mid_disp','spec','bunching'],s_inp=['max'],z_inp=[0,'end'], savefig=1, saveval=1, show=1)
t4 =time.time()
delta2 = t4-t3
print 'total post-processing time = ',delta2